function [a_sims_all,chat_sims_all,cstar_sims_all,aphat_sims_all,astar_sims_all,eta_sims_all,sigma_etasq_sims_all,alpha_star_sims_all,mpc_all,mpc_all2, mpc_all3, mpc_all4, age_sims_all,prior_sims_all, chat_prime_all, alpha_mpc_sims_all]...
    = bc_endsignals_sim(w,param,y_sims_all,eps_eta_mat,reset_shocks,b,mu0_fun, ctilde, stateGrid,stateGrid_step,stateGrid_size,cstar, log_c, a_sims0, ex_signals)

% INPUT
%   param:          parameter vector
%                   [bta,r,abar,W_cc,kappa,delta,deltax,ybar,sigma_y,
%                   cstar_ya,cbar,prior_mean_coeff,sigma_c,psi,cov_fun]
%   y_sims_all:     simulated y                             [sim_periods,nsim]
%   eps_eta_mat:    eta innovations (epsilon_eta)           [sim_periods,nsim]
%   reset_shocks:   reset shocks                            [sim_periods,nsim]
%   b:              borrowing constraint (<=0)              [1,1]
%   tilde           whether to use ctilde (T) or cstar (F)  boolean

% OUTPUT
%   a_sims_all:             simulated values for a              [sim_periods,nsim]
%   chat_sims_all:          simulated values for chat           [sim_periods,nsim]
%   cstar_sims_all:         simulated values for cstar          [sim_periods,nsim]
%   aphat_sims_all:         simulated values for aphat          [sim_periods,nsim]
%   apstar_sims_all:        simulated values for apstar         [sim_periods,nsim]
%   eta_sims_all:           simultaed values for eta            [sim_periods,nsim]
%   sigma_etasq_sims_all:   simulated values for sigma_etasq    [sim_periods,nsim]

% unpack parameters
bta = param(1);
r = param(2);
abar = param(3);
W_cc = param(4);
kappa = param(5);
deltax = param(6);
sigma_c = param(7);
psi = param(8);
cov_fun = param(9);
gam = param(10);
sigma_l = param(11);
lbar = param(12);
delta_eta = param(13);
delta_inh= param(14);
factor_mpc1 = param(15);
factor_mpc2 = param(16);
delta_eta_infoshock = param(17);
psi_tau  = 0;     

[sim_periods, nsim] = size(y_sims_all);



a_sims_all          = NaN(sim_periods,nsim);
chat_sims_all       = NaN(sim_periods,nsim);
cstar_sims_all      = NaN(sim_periods,nsim);
aphat_sims_all      = NaN(sim_periods,nsim);
astar_sims_all     = NaN(sim_periods,nsim);
eta_sims_all        = NaN(sim_periods,nsim);
sigma_etasq_sims_all= NaN(sim_periods,nsim);
chat_prime_all      = NaN(sim_periods,nsim);

alpha_star_sims_all  = NaN(sim_periods,nsim);
alpha_mpc_sims_all   = NaN(sim_periods, nsim);
mpc_all              = NaN(sim_periods,nsim);
mpc_all2             = NaN(sim_periods,nsim);
mpc_all3             = NaN(sim_periods,nsim);
mpc_all4             = NaN(sim_periods,nsim);

age_sims_all        = NaN(sim_periods,nsim);
prior_sims_all      = NaN(sim_periods,nsim);

parfor j=1:nsim
    
    psi_tau = 0; 

    %  preallocate
    y_sims           = y_sims_all(:,j);
    a_sims           = NaN(sim_periods,1);
    chat_sims        = NaN(sim_periods,1);
    cstar_sims       = NaN(sim_periods,1);
    aphat_sims       = NaN(sim_periods,1);
    astar_sims       = NaN(sim_periods,1);
    eta_sims         = NaN(sim_periods,1);
    sigma_etasq_sims = NaN(sim_periods,1);
    sigma_etasq_save = NaN(sim_periods,1);
    eps_eta_sims     = eps_eta_mat(:,j);
    chat_prime_sims  = NaN(sim_periods,1);

    alpha_star_sims = NaN(sim_periods,1);
    alphat_mpc1_sims = NaN(sim_periods,1);
    age_sims     = NaN(sim_periods,1);
    prior_sims     = NaN(sim_periods,1);
    mpc_sims              =NaN(sim_periods,1);
    mpc_sims2              =NaN(sim_periods,1);
    mpc_sims3              =NaN(sim_periods,1);
    mpc_sims4              =NaN(sim_periods,1);

    %  initial states
    a_sims(1) = a_sims0(j);
    mu0_ya_hist = [];
    mu0_ya_hist_mpc = [];
    hist_mat = [];
    Linv = [];
    mu_t0_hist = [];


    %  period 1
    [alpha_t, prior_mean ,chat, eta_t, sigma_etasq, mu0_ya_hist, hist_mat, Linv, mu_t0_hist] = bc_update_iid_fast(r,  W_cc, kappa, w*y_sims(1),         a_sims(1),         sigma_c, psi, cov_fun, ex_signals, eps_eta_sims(1),         b, mu0_ya_hist, hist_mat,Linv,mu_t0_hist, ctilde, stateGrid,stateGrid_step,stateGrid_size, log_c, psi_tau);
    if log_c == 1
        chat = exp(chat);
    end
    chat_sims(1) = chat;
    aphat_sims(1) = w*y_sims(1)+ (1+r)*a_sims(1)-chat_sims(1);

    eta_sims(1)  = eta_t;
    sigma_etasq_sims(1) = sigma_etasq;
    sigma_etasq_save(1) = sigma_etasq;

    alpha_star_sims(1)=alpha_t;
    prior_sims(1)=prior_mean;
    age_sims(1)=1;

    cstar_sims(1)    = min(cstar(a_sims(1)*(1+r) + w*y_sims(1)), a_sims(1)*(1+r) + w*y_sims(1) - b);

    prior_prime = bc_update_mut0(r, W_cc, kappa, w*y_sims(2), w*y_sims(1) + (1+r)*a_sims(1) - chat_sims(1), sigma_c, psi, cov_fun, 0, eps_eta_sims(1),  b, mu0_ya_hist, hist_mat,Linv,mu_t0_hist, ctilde, stateGrid,stateGrid_step,stateGrid_size, log_c);
    chat_prime_sims(2) = prior_prime;


    %i= 1;
    y_mpc= y_sims(1) + factor_mpc1*sigma_l;
    [~,~,chat_mpc, eta_temp, ~, ~, ~, ~, ~] = bc_update_iid_fast(r, W_cc, kappa, w*y_mpc, a_sims(1), sigma_c, psi, cov_fun, ex_signals, eps_eta_sims(1),  b, mu0_ya_hist, hist_mat,Linv,mu_t0_hist, ctilde, stateGrid,stateGrid_step,stateGrid_size, log_c,psi_tau);
    
    mpc_sims(1)=(chat_mpc-chat_sims(1))/(w*factor_mpc1*sigma_l);

    if log_c == 1
        mpc_sims(1)=(exp(chat_mpc)-chat_sims(1))/(w*factor_mpc1);
    end

    %%%% big positive shock
    y_mpc= y_sims(1) + factor_mpc2*sigma_l;
    [~,~,chat_mpc, ~, ~, ~, ~, ~, ~] = bc_update_iid_fast(r, W_cc, kappa, w*y_mpc, a_sims(1), sigma_c, psi, cov_fun, ex_signals, 0*eps_eta_sims(1),  b, mu0_ya_hist, hist_mat,Linv,mu_t0_hist, ctilde, stateGrid,stateGrid_step,stateGrid_size, log_c,psi_tau);
    mpc_sims2(1)=(chat_mpc-chat_sims(1))/(w*factor_mpc1);
    if log_c == 1
        mpc_sims2(1)=(exp(chat_mpc)-chat_sims(1))/(w*factor_mpc1);
    end

    

    %  period 2+
    for i = 2:sim_periods

        %if i == 1301
        %    break
        %end

        if reset_shocks(i, j) < deltax

            

            sigma_etasq_sims(1:i-1) = sigma_etasq_sims(1:i-1);
            a_sims(i)     = (w*y_sims(i-1) + (1+r)*a_sims(i-1) - chat_sims(i-1))*delta_inh;
            %astar_sims(i) = w*y_sims(i-1) + (1+r)*a_sims(i-1) - cstar_sims(i-1);
            astar_sims(i) = a_sims(i);
            %%%% Carry over beliefs properly


            if sum(reset_shocks(i,:) < deltax ,2) == nsim
                %sigma_etasq_sims(1:i-1) = sigma_etasq_sims(1:i-1)/delta_eta_infoshock;
                sigma_etasq_sims(1:i-1) = sigma_etasq_sims(1:i-1);

            else

                sigma_etasq_sims(1:i-1) = sigma_etasq_sims(1:i-1)/delta_eta;

            end

            sigma_etasq_hist = sigma_etasq_sims( sigma_etasq_sims < 1e10);
            eta_hist         = eta_sims( sigma_etasq_sims < 1e10);
            ya_hist          = y_sims( sigma_etasq_sims < 1e10) + (1+r)*a_sims( sigma_etasq_sims < 1e10);

            [mu_t0_hist, ~] = gp_update_tilde(ya_hist, ya_hist, eta_hist, sigma_etasq_hist, mu0_fun,sigma_c,psi, cov_fun);


            mu0_ya_hist = mu0_fun( ya_hist );


            hist_mat = [ y_sims( sigma_etasq_sims < 1e10), a_sims( sigma_etasq_sims < 1e10), eta_hist, sigma_etasq_hist];

            temp1 = ya_hist - ya_hist';
            Sigma_y1yN = sigma_c^2*exp(-psi*(temp1.^2));
            Sigma22 = Sigma_y1yN + diag(ones(length(ya_hist),1).*(sigma_etasq_hist));

            L = chol(Sigma22, 'lower');

            Linv = L\eye( length(Sigma22));
            %}

            age_sims(i)=1;

        else
            %  discount past info: incr variance of eta (decr precision of eta)
            sigma_etasq_sims(1:i-1) = sigma_etasq_sims(1:i-1);
            %  a_t
            a_sims(i)    = w*y_sims(i-1) + (1+r)*a_sims(i-1) - chat_sims(i-1);
            astar_sims(i) = w*y_sims(i-1) + (1+r)*a_sims(i-1) - cstar_sims(i-1);
            age_sims(i)=age_sims(i-1)+1;


        end


        if sum(reset_shocks(i,:) < deltax ,2) == nsim
            psi_tau = delta_eta_infoshock;
        end

        %small positive shock
        y_mpc= y_sims(i) + factor_mpc1*sigma_l;
        [alphat_mpc1,~,chat_mpc1, ~, ~, ~, ~, ~, ~] = bc_update_iid_fast(r, W_cc, kappa, w*y_mpc, a_sims(i), sigma_c, psi, cov_fun, 1e3*ex_signals, 1*eps_eta_sims(i),  b, mu0_ya_hist, hist_mat,Linv,mu_t0_hist, ctilde, stateGrid,stateGrid_step,stateGrid_size, log_c,psi_tau);



        %%%% y/12 negative
        %large positive shock
        y_mpc = y_sims(i) + factor_mpc2*sigma_l;
        [~,~,chat_mpc2, ~, ~, ~, ~, ~, ~] = bc_update_iid_fast(r, W_cc, kappa, w*y_mpc, a_sims(i), sigma_c, psi, cov_fun, 1e3*ex_signals, 1*eps_eta_sims(i),  b, mu0_ya_hist, hist_mat,Linv,mu_t0_hist, ctilde, stateGrid,stateGrid_step,stateGrid_size, log_c,psi_tau);



        %  chat_t
        [alpha_t,prior_mean,chat, eta_t, sigma_etasq, mu0_ya_hist, hist_mat,Linv,mu_t0_hist] = bc_update_iid_fast(r, W_cc, kappa, w*y_sims(i), a_sims(i), sigma_c, psi, cov_fun, ex_signals, eps_eta_sims(i),  b, mu0_ya_hist, hist_mat,Linv,mu_t0_hist, ctilde, stateGrid,stateGrid_step,stateGrid_size, log_c,psi_tau);
        if log_c == 1
            chat = exp(chat);
        end
        chat_sims(i)     = chat;
        cstar_sims(i)    = min(cstar(astar_sims(i)*(1+r) + w*y_sims(i)), astar_sims(i)*(1+r) + w*y_sims(i) - b);
        aphat_sims(i)    = w*y_sims(i)+ (1+r)*a_sims(i)-chat_sims(i);

        if i < sim_periods
            prior_prime  = bc_update_mut0(r, W_cc, kappa, w*y_sims(i+1), w*y_sims(i) + (1+r)*a_sims(i) - chat_sims(i), sigma_c, psi, cov_fun, 0, eps_eta_sims(i),  b, mu0_ya_hist, hist_mat,Linv,mu_t0_hist, ctilde, stateGrid,stateGrid_step,stateGrid_size, log_c);
            chat_prime_sims(i+1) = prior_prime;
        end

        eta_sims(i)  = eta_t;  % signal at t
        sigma_etasq_sims(i) = sigma_etasq;  % precision of signal at t
        sigma_etasq_save(i) = sigma_etasq;

        alpha_star_sims(i)=alpha_t;
        prior_sims(i)=prior_mean;
        alphat_mpc1_sims(i) = alphat_mpc1;


        
        %%%% y/12 positive
        mpc_sims(i)=(chat_mpc1-chat_sims(i))/(w*factor_mpc1*sigma_l);
        if log_c == 1
            mpc_sims(i)=(exp(chat_mpc1)-chat_sims(i))/(w*factor_mpc1*y_sims(i));
        end

        mpc_sims2(i)=(chat_mpc2-chat_sims(i))/(w*factor_mpc2*sigma_l);
        if log_c == 1
            mpc_sims2(i)=(exp(chat_mpc2)-chat_sims(i))/(w*factor_mpc2*y_sims(i));
        end


        %%%% y/4 positive
        y_mpc=y_sims(i) + factor_mpc1*sigma_l;
        [a_temp2,~,chat_mpc3, ~, ~, ~, ~, ~, ~] = bc_update_iid_fast(r, W_cc, 1000*kappa, w*y_mpc, a_sims(i), sigma_c, psi, cov_fun, 1e3*ex_signals, 0*eps_eta_sims(i),  b, mu0_ya_hist, hist_mat,Linv,mu_t0_hist, ctilde, stateGrid,stateGrid_step,stateGrid_size, log_c,psi_tau);


        %%%% y/4 negative
        y_mpc=y_sims(i) +  factor_mpc2*sigma_l;
        [~,~,chat_mpc4, ~, ~, ~, ~, ~, ~] = bc_update_iid_fast(r, W_cc, kappa, w*y_mpc, a_sims(i), sigma_c, psi, cov_fun, 1e3*ex_signals, 0*eps_eta_sims(i),  b, mu0_ya_hist, hist_mat,Linv,mu_t0_hist, ctilde, stateGrid,stateGrid_step,stateGrid_size, log_c,psi_tau);


        mpc_sims3(i)=(chat_mpc3-chat_sims(i))/(w*factor_mpc1*sigma_l);
        if log_c == 1
            mpc_sims3(i)=(exp(chat_mpc3)-chat_sims(i))/(w*factor_mpc1*y_sims(i));
        end


        mpc_sims4(i)=(chat_mpc4-chat_sims(i))/(w*factor_mpc2*sigma_l);
        if log_c == 1
            mpc_sims4(i)=(exp(chat_mpc4)-chat_sims(i))/(w*factor_mpc2*y_sims(i));
        end


    end

    a_sims_all(:,j)         = a_sims;
    chat_sims_all(:,j)      = chat_sims;
    cstar_sims_all(:,j)     = cstar_sims;
    aphat_sims_all(:,j)     = aphat_sims;
    astar_sims_all(:,j)     = astar_sims;
    alpha_mpc_sims_all(:,j)    = alphat_mpc1_sims;
    eta_sims_all(:,j)       = eta_sims;
    sigma_etasq_sims_all(:,j) = sigma_etasq_save;
    chat_prime_all(:,j)   = chat_prime_sims;

    alpha_star_sims_all(:,j)= alpha_star_sims;
    age_sims_all(:,j)= age_sims;
    prior_sims_all(:,j)= prior_sims;
    mpc_all(:,j)=mpc_sims;
    mpc_all2(:,j)=mpc_sims2;
    mpc_all3(:,j)=mpc_sims3;
    mpc_all4(:,j)=mpc_sims4;

end

end